Phase Diagram of Interacting Bosons on the Honeycomb Lattice 
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We study the ground state properties of repulsively interacting bosons on the honeycomb lattice 
using large-scale quantum Monte Carlo simulations. In the hard-core limit the half-filled system 
develops long ranged diagonal order for sufficiently strong nearest-neighbor repulsion. This staggered 
solid melts at a first order quantum phase transition into the superfluid phase, without the presence 
of any intermediate supersolid phase. Within the superfluid phase, both the superfluid density and 
the compressibility exhibit local minima near particle- (hole-) density one quarter, while the density 
and the condensate fraction show inflection points in this region. Relaxing the hard-core constraint, 
supersolid phases emerge for soft-core bosons. The suppression of the superfluid density is found to 
persist for sufficiently large, finite on-site repulsion. 

PACS numbers: 05.30.Jp, 03.75.Lm, 75.10.Jm, 75.40.Mg 
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The properties of interacting bosons confined to peri- 
odic lattice structures are intensively being studied re- 
cently. This effort is driven by both the progress in re- 
alizing such many-body systems using ultra-cold atoms 
on optical lattices [l], d, y|, 3, as well as the interest- 
ing phases and quantum phase transitions expected to 
emerge. 

In particular, the ground state can be superfluid, 
Mott insulating or, in the presence of disorder, a Bose 
glass 0, 0, 0- Recently, based on unbiased numerical 
simulations, also supersolid [1, & EHE, 12, 13, 14, 15 1 
and exotic valence-bond-solid phases [1 61 ] were confirmed 
to be acessible, and unconventional quantum critical- 
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was discussed for such systems, 



in particular on non-bipartite, geometrically frustrated 
two-dimensional lattices, such as the triangular or the 
Kagomc lattice. 

Concerning bosons on two-dimensional bipartite, and 
thus non-frustrated lattices, so far mainly a square lat- 
tice structure has also been considered in various other 
studies [19-31]. However, in light of the recent real- 
ization of two-dimensional graphene, with its noticeable 
properties of fermions on the underlying honeycomb lat- 
tice [32l. l33l. l34l. l35j . the question arrises, how interacting 
bosons behave on the honeycomb lattice. While both the 
honeycomb and the square lattice are bipartite, quantum 
flucuations are expected to be more relevant on the hon- 
eycomb lattice, due to its lowest possible coordination in 
two dimensions. 

Motivated by these considerations, we analyze in this 
Letter the ground state properties of the extended Bose- 
Hubbard model 
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on the honeycomb lattice. Here, b\ (pi) denote bosonic 



FIG. 1: (Color online) Ground state phase diagram of bosons 
on the honeycomb lattice in the hard-core limit. The black 
dot locates the Heisenberg point, where the model has an 
enhanced SU (2) symmetry. In the inset the superfluid density 
ps and the structure factor Sa are shown along t/V = 0.1 
(indicated by the vertical bar in the main part of the figure) . 



creation (annihilation) operators for bosons on lattice 
site i, t is the nearest-neighbor hopping amplitude, U an 
onsite repulsion, V a nearest-neighbor repulsion, and (i 
the chemical potential in the grand-canonical ensemble, 
which controls the filling of the lattice. 

The honeycomb lattice is bipartite, with a two-site unit 
cell and a uniform coordination number z = 3; the unit 
cell contains one site from each sublattice. In the fol- 
lowing, we employ the stochastic series expansion 36 1 
quantum Monte Carlo (QMC) method with directed loop 
updates [37], [3], and study the Hamiltonian in Eq. (1) 
on finite lattices with N = 2L 2 sites for L up to 32, at 
temperatures sufficiently low in order to resolve ground 
state properties of these finite systems § . In contrast to 
geometrically frustrated lattices, the sign of t in Eq. (1) 
is irrelevant on the bipartite honeycomb lattice, in that 
also t < would cause no QMC sign problem. 
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FIG. 2: (Color online) Density p (with respect to half-filling) 
and condensate fraction po/p at t/V — 0.1 as functions of 
(i/V for bosons on the honeycomb lattice in the hard-core 
limit. The inset shows the compressibility k which exhibits a 
local minimum in the superfluid regime. 



We first consider the hard-core limit, U/t — > oo, of 
Eq. (1), for which the bosonic model can be mapped 
onto the spin- 1/2 XXZ model (3^|. This allows for an 
interpretation of our results in terms of both bosons and 
quantum spins. The ground state phase diagram of the 
grand-canonical ensemble in the hard-core limit is shown 
in Fig. 1. Considering the one-particle and one-hole prob- 
lem, one finds that for p < —zt the system is empty, and 
for p/V > z(l + t/V) it is fully occupied by one bo- 
son per lattice site (density p — 1). At large values of 
t/V > 0.5, the bosons form a superfluid (SF) with off- 
diagonal long-range order (ODLRO), characterized by a 
finite superfluid density, which in the QMC simulations 
is obtained as p s = (W 2 ) / (2/3t) from the particle winding 
number fluctuations (W 2 ), where (3 denotes the inverse 
temperature [Z(|. Furthermore, for sufficiently small val- 
ues of t/V, the system exhibits a solid phase of density 
p = 1/2. Its largest extend is set by the Hcisenberg 
point (t/V,/j,/V) = (1/2,3/2), where the Hamiltonian in 
Eq. (1) has an enhanced SU(2) symmetry. In the perfect 
p = 1/2 solid at t = 0, one out of the two sublattices on 
the honeycomb lattice is occupied, the other being empty 
(so that reflection symmetry, but not lattice translational 
symmetry is broken) . The corresponding structure factor 
for this alternating (staggered) diagonal long-range order 
(DLRO) is thus given by 



(2) 



where e\ = ±1 for i on sublattice A(B), in terms of the 
density-density correlations. The inset of Fig. 1 shows 
the evolution of both ps and the order parameter Sa/N 
along a cut through the phase diagram at t/V = 0.1 (the 
vertical bar in the main part of Fig. 1), and corresponds 
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FIG. 3: (Color online) Superfluid density of bosons on the 
honeycomb lattice at t/V = 0.2 as functions of p/V for vari- 
ous values of U/t, and in the hard-core limit. 



to their thermodynamic limit behavior. Both quantities 
clearly exhibit the first-order nature of the solid-SF 
quantum phase transition, and the absence of any 
intermediate supcrsolid (SS) phase with simultaneous 
DLRO and ODLRO (only at the Hcisenberg point 
may DLRO and ODLRO co-exist, due to the SU(2) 
symmetry). In the hard-core limit, the phase diagram 
of Eq. (1) is thus similar to the one of hard-core bosons 
on the square lattice [1], irrespectively a rescaled /i-axis, 
due to the different coordination number (z = 4 on the 
square lattice). In both cases, domain wall proliferation 
at the quantum melting transition of the half-filled solid 
renders a SS unstable towards phase separation in the 
hard-core limit of Eq. (1) Note, that also through the 
Heisenberg point the quantum melting of the p = 1/2 
solid is discontinuous [29L l3ll . 41 1, despite attempts to 
extract critical behavior from finite size scaling plots of 
PS [Hj E3l- In fact) the extracted effective dynamical 
critical exponent z decreases to zero, taking sufficiently 
large system sizes [4l|, concistent with the trend ob- 
served by comparing the values obtained thus far [2^, \^ . 



In contrast to the case of the square lattice, we observe 
on the honeycomb lattice a pronounced dip of pg inside 
the SF phase (c.f. the inset of Fig. 1), with a local 
minimum at a filling pu close to 3/4, corresponding to 
a hole density of 1/4 (p D w 0.744 at (n/V) D » 2.97 
for t/V = 0.1). Due to particle- hole symmetry in the 
hard-core limit, this suppression of ps also appears near 
a particle filling of p = 1/4. Fig. 3 shows a similar dip 
in ps for t/V = 0.2, where we again find po close to 
3/4. This indicates the presence of geometric hindrance 
in the superfluid flow on the honeycomb lattice, well 
inside the SF region. Indeed, also the compressibility 
k = dp/dp: exhibits a minimum at (p/V)o, c.f. the 
inset of Fig. 2. In contrast, as seen in Fig. 2, both 
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FIG. 4: (Color online) Ground state phase diagram of bosons 
on the honeycomb lattice for finite U/t = 20, with supersolid 
(SS) phases emerging upon doping the solid phases. Solid 
(dashed) lines denote continuous (first-order) quantum phase 
transitions. 
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FIG. 5: (Color online) Density p, superfluid density p a , and 
structure factor Sa as functions of p/V for bosons on the 
honeycomb lattice for U/t = 20 along t/V = 0.14 (upper 
panel) and t/V = 0.16 (lower panel). 



the density p and the condensate fraction, po/p, where 
Po = (b ' (k = 0)6(k = 0)), increase in that region upon 
increasing p, and show an inflection point at (p/V)o- 
We verified that near the dip the system does not 
exhibit any incommensurate long-range order neither 
in the density-density nor in bond-bond correlation 
functions [lit . 

Extending the analysis beyond the hard-core limit, we 
next consider soft-core bosons on the honeycomb lattice, 
and assess the presence of SS phases in this model. Fig. 4 
shows the ground state phase diagram for U/t = 20 inside 
the regime where p < 1. The phase diagram exhibits sev- 
eral additional phases; in particular, a SS phase emerges 
upon doping the p = 1/2 solid with additional bosons. 
The simultaneous presence of DLRO and ODLRO within 
the SS regime for t/V = 0.14 and t/V = 0.16 can be seen 
from Fig. 5, by finite values of both Sa/N and ps inside 
a range of p/V values. For (zV — U)~t, the additional 
particles form a superfluid atop the p = 1/2 solid back- 
ground, due to a kinetic energy gain which prevents the 
domain wall proliferation Hole doping of the p = 1/2 
solid however does not result in a similar kinetic energy 
gain, and thus no SS state emerges for p < 1/2. 

The nature of the incompressible state with p = 1 de- 
pends on the ratio zV/U [8(: In the soft-core case, taking 
U/t = 20, we find for t/V < 0.15 (i.e., zV/U > 1), that 
the system forms a p = 1 solid with a finite value of Sa/N 
in the thermodynamic limit. This p = 1 solid corresponds 
to each site of one out of the two sublattices being oc- 
cupied by two bosons, the other sublattice being empty. 
In contrast, for t/V > 0.15 (i.e., zV/U < 1), the sys- 
tem is a uniform p = 1 Mott insulator, with Sa/N = 0. 
This different behavior at p = 1 is seen for t/V = 0.14 



and t/V = 0.16 in Fig. 5, and results from the com- 
petition between onsite- and nearest-neighbor repulsion 
terms in Eq. (1) Fig. 5 furthermore indicates, that 
the transition from the p < 1 SS to the p = 1 solid 
is strongly first-order, whereas the SS-SF transition and 
the SF-Mott transition are continuous, as expected also 
from kinetic energy considerations Q. As indicated in 
Fig. 4, a further SS phase, with p > 1, emerges out of 
the p = 1 solid phase upon further increasing p. Details 
of the phase diagram for fillings p > 1 will be presented 
elsewhere (43| . 

For U/t = 20, the superfluid density curves in Fig. 5 
do not exhibit local minima such as the pronounced dip 
in ps observed in the hard-core limit (in the inset of 
Fig. 1). In order to assess, whether this suppression is 
restricted to the hard-core limit of Eq. (1), we performed 
simulations also for increasing values of U/t — 30,50, 
and 80. In Fig. 3, the resulting values of ps are shown 
as functions of p/V for a common value of t/V = 0.2, 
and compared to the hard-core limit. We find that 
for U/t = 30 the superfluid density is almost flat near 
p/V w 2.9, and for sufficiently large onsite repulsion U, 
it indeed develops a shallow local minimum, similar to 
the dip in the hard-core limit. The reduced superfluidity 
thus relates to the interplay of a large energy penalty 
for double occupation, and the low connectivity on the 
honeycomb lattice, compared to the square lattice. In 
fact, we find that on (inhomogencous) two-dimensional 
lattices with an even lower average coordination number 
z < 3, the superfluid density gets reduced to zero, and 
incompressible phases with fractional fillings p = 1/4 
and 3/4 emerge [43( | . 

In conclusion, we studied the ground state properties 
of bosons on the honeycomb lattice, as described by 
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the extended Bosc-Hubbard model. We found that 
the phases and the structure of the phase diagram 
correspond to those of the same model on the square 
lattice, both in the hard-core limit and for finite on-site 
repulsion U. This is anticipated given the bipartiteness 
of both lattice structures. In particular, a staggered 
solid phase exists at half-filling for sufficiently strong 
nearest-neighbor repulsion V. Doping this solid with 
additional bosons leads to a supersolid phase only 
for finite values of U, whereas in the hard-core limit 
and for hole-doping the system phase separates. The 
incompressible state at filling one is a Mott insulator for 
dominant U > zV, and a staggered solid for U < zV. 
Furthermore, we found that for sufficiently large values 
of U, both the superfluid density and the compressibility 
are suppressed inside the superfluid phase at particle 
fillings of about one quarter and three quarters. As the 
condensate density does not show a similar reduction, 
bosons on the honeycomb lattice thus feature regions, 
where an increase of the condensate in the bosonic 
system contrasts to a decrease in its superfluid response. 

We would like to thank F. Alct, R. Moessner, and A. 
Muramatsu for helpful discussions, and NIC Jiilich and 
HRLS Stuttgart for allocation of CPU time. 

Note added.- After finishing the numerical calculations, 
we became aware of a parallel work |42j |. where results 
partially similar to our findings are presented. 
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